home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / msdos / raytrace / rayshade / src / cone.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-07-18  |  6.9 KB  |  321 lines
  1. /*
  2.  * cone.c
  3.  *
  4.  * Copyright (C) 1989, 1991, Craig E. Kolb
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * $Id: cone.c,v 4.0 91/07/17 14:36:46 kolb Exp Locker: kolb $
  17.  *
  18.  * $Log:    cone.c,v $
  19.  * Revision 4.0  91/07/17  14:36:46  kolb
  20.  * Initial version.
  21.  * 
  22.  */
  23. #include "geom.h"
  24. #include "cone.h"
  25.  
  26. static Methods *iConeMethods = NULL;
  27. static char coneName[] = "cone";
  28.  
  29. unsigned long ConeTests, ConeHits;
  30.  
  31. Cone *
  32. ConeCreate(br, bot, ar, apex)
  33. Vector *bot, *apex;
  34. Float br, ar;
  35. {
  36.     Cone *cone;
  37.     Float tantheta, lprime, tlen, len, dtmp;
  38.     Vector axis, base, tmp;
  39.  
  40.     /*
  41.      * The passed basepoint must be closer to the origin of the
  42.      * cone than the apex point, implying that the base radius
  43.      * must be smaller than the apex radius.  If the values passed
  44.      * reflect the opposite, we switch everything.
  45.      */
  46.     if(ar < br) {
  47.         tmp = *bot;
  48.         *bot = *apex;
  49.         *apex = tmp;
  50.         dtmp = br;
  51.         br = ar;
  52.         ar = dtmp;
  53.     } else if (equal(ar, br)) {
  54.         RLerror(RL_WARN, "Cone is a cylinder -- discarding.\n");
  55.         return (Cone *)NULL;
  56.     }
  57.     /*
  58.      * Find the axis and axis length.
  59.      */
  60.     VecSub(*apex, *bot, &axis);
  61.     len = VecNormalize(&axis);
  62.     if (len < EPSILON) {
  63.         RLerror(RL_WARN, "Degenerate cone.\n");
  64.         return (Cone *)NULL;
  65.     }
  66.  
  67.     cone = (Cone *)share_malloc(sizeof(Cone));
  68.  
  69.     /*
  70.      * To render a cone, we transform the desired cone into
  71.      * a canonical, Z-axis aligned, unit length, unit radius
  72.      * at the apex cone.
  73.      *
  74.      * Here, we construct the transformation matrix to map
  75.      * from cone<-->world space.
  76.      */
  77.  
  78.     /*
  79.      * "tantheta" is the change in radius per unit length along
  80.      * the cone axis.
  81.      */
  82.     tantheta = (ar - br) / len;
  83.     /*
  84.      * lprime defines the distance along the axis where the cone starts
  85.      */
  86.     lprime = br / tantheta;
  87.     /*
  88.      * Find the true base (origin) of the cone.
  89.      */
  90.     VecScale(-lprime, axis, &base);
  91.     VecAdd(base, *bot, &base);
  92.     /*
  93.      * tlen is the total length of the cone.
  94.      */
  95.     tlen = lprime + len;
  96.     /*
  97.      * start_dist is the distance from the origin of the canonical
  98.      * cone at which the cone begins
  99.      */
  100.     cone->start_dist = lprime / tlen;
  101.     CoordSysTransform(&base, &axis, ar, tlen, &cone->trans);
  102.     return cone;
  103. }
  104.  
  105. Methods *
  106. ConeMethods()
  107. {
  108.     if (iConeMethods == (Methods *)NULL) {
  109.         iConeMethods = MethodsCreate();
  110.         iConeMethods->name = ConeName;
  111.         iConeMethods->create = (GeomCreateFunc *)ConeCreate;
  112.         iConeMethods->methods = ConeMethods;
  113.         iConeMethods->intersect = ConeIntersect;
  114.         iConeMethods->normal = ConeNormal;
  115.         iConeMethods->uv = ConeUV;
  116.         iConeMethods->bounds = ConeBounds;
  117.         iConeMethods->stats = ConeStats;
  118.         iConeMethods->checkbounds = TRUE;
  119.         iConeMethods->closed = FALSE;
  120.     }
  121.     return iConeMethods;
  122. }
  123.  
  124. /*
  125.  * Ray-cone intersection test.  This routine is far from optimal, but
  126.  * it's straight-forward and it works...
  127.  */
  128. int
  129. ConeIntersect(cone, ray, mindist, maxdist)
  130. Cone *cone;
  131. Ray *ray;
  132. Float mindist, *maxdist;
  133. {
  134.     Float t1, t2, a, b, c, disc, zpos, distfact;
  135.     Ray newray;
  136.     Vector nray, npos;
  137.     Float nmin;
  138.  
  139.     ConeTests++;
  140.  
  141.     /*
  142.      * Transform ray from world to cone space.
  143.      */
  144.     newray = *ray;
  145.     distfact = RayTransform(&newray, &cone->trans.itrans);
  146.     nray = newray.dir;
  147.     npos = newray.pos;
  148.     nmin = mindist * distfact;
  149.  
  150.     a = nray.x * nray.x + nray.y * nray.y - nray.z*nray.z;
  151.     b = nray.x * npos.x + nray.y * npos.y - nray.z*npos.z;
  152.     c = npos.x*npos.x + npos.y*npos.y - npos.z*npos.z;
  153.  
  154.     if (fabs(a) < EPSILON) {
  155.         /*
  156.          * Only one intersection point...
  157.          */
  158.         t1 = -0.5*c / b;
  159.         zpos = npos.z + t1 * nray.z;
  160.         if (t1 < nmin || zpos < cone->start_dist || zpos > 1.)
  161.             return FALSE;
  162.         t1 /= distfact;
  163.         if (t1 < *maxdist) {
  164.             *maxdist = t1;
  165.             ConeHits++;
  166.             return TRUE;
  167.         }
  168.         return FALSE;
  169.     } else {
  170.         disc = b*b - a*c;
  171.         if (disc < 0.)
  172.             return FALSE;        /* No possible intersection */
  173.         disc = sqrt(disc);
  174.         t1 = (-b + disc) / a;
  175.         t2 = (-b - disc) / a;
  176.         /*
  177.          * Clip intersection points.
  178.          */
  179.         zpos = npos.z + t1 * nray.z;
  180.         if (t1 < nmin || zpos < cone->start_dist || zpos > 1.) {
  181.             zpos = npos.z + t2 * nray.z;
  182.             if (t2 < nmin || zpos < cone->start_dist ||
  183.                 zpos > 1.)
  184.                 return FALSE;
  185.             else
  186.                 t1 = t2 / distfact;
  187.         } else {
  188.             zpos = npos.z + t2 * nray.z;
  189.             if (t2 < nmin || zpos < cone->start_dist ||
  190.                 zpos > 1.)
  191.                 t1 /= distfact;
  192.             else
  193.                 t1 = min(t1, t2) / distfact;
  194.         }
  195.         if (t1 < *maxdist) {
  196.             *maxdist = t1;
  197.             ConeHits++;
  198.             return TRUE;
  199.         }
  200.         return FALSE;
  201.     }
  202. }
  203.  
  204. /*
  205.  * Compute the normal to a cone at a given location on its surface.
  206.  */
  207. int
  208. ConeNormal(cone, pos, nrm, gnrm)
  209. Cone *cone;
  210. Vector *pos, *nrm, *gnrm;
  211. {
  212.     Vector npos;
  213.  
  214.     /*
  215.      * Transform intersection point to cone space.
  216.      */
  217.     npos = *pos;
  218.     PointTransform(&npos, &cone->trans.itrans);
  219.     
  220.     /*
  221.      * The following is equal to
  222.      * (pos X (0, 0, 1)) X pos
  223.      */
  224.     nrm->x = npos.x * npos.z;
  225.     nrm->y = npos.y * npos.z;
  226.     nrm->z = -npos.x * npos.x - npos.y * npos.y;
  227.  
  228.     /*
  229.      * Transform normal back to world space.
  230.      */
  231.     NormalTransform(nrm, &cone->trans.itrans);
  232.     *gnrm = *nrm;
  233.     return FALSE;
  234. }
  235.  
  236. void
  237. ConeUV(cone, pos, norm, uv, dpdu, dpdv)
  238. Cone *cone;
  239. Vector *pos, *norm, *dpdu, *dpdv;
  240. Vec2d *uv;
  241. {
  242.     Vector npos;
  243.     Float val;
  244.  
  245.     npos = *pos;
  246.     PointTransform(&npos, &cone->trans.itrans);
  247.  
  248.     uv->v = (npos.z - cone->start_dist) / (1. - cone->start_dist);
  249.     if (npos.z < EPSILON)
  250.         uv->u = 0.;
  251.     else {
  252.         val = npos.x / npos.z;
  253.         /*
  254.          * Be careful not to feed |val| > 1 to acos
  255.          */
  256.         if (val > 1.)
  257.             uv->u = 0.;
  258.         else if (val < -1.)
  259.             uv->u = 0.5;
  260.         else {
  261.             uv->u = acos(val) / TWOPI;
  262.             if (npos.y < 0.)
  263.                 uv->u = 1. - uv->u;
  264.         }
  265.     }
  266.     /*
  267.      * dpdv = pos
  268.      * dpdu = dpdv X norm = (-norm->y, norm->x, 0.)
  269.      */
  270.     if (dpdu) {
  271.         *dpdv = npos;
  272.         dpdu->x = -npos.y;
  273.         dpdu->y = npos.x;
  274.         dpdu->z = 0.;
  275.         VecTransform(dpdu, &cone->trans.trans);
  276.         VecTransform(dpdu, &cone->trans.trans);
  277.         (void)VecNormalize(dpdu);
  278.         (void)VecNormalize(dpdv);
  279.     }
  280. }
  281.  
  282. /*
  283.  * Return the extent of a cone.
  284.  */
  285. void
  286. ConeBounds(cone, bounds)
  287. Cone *cone;
  288. Float bounds[2][3];
  289. {
  290.     bounds[LOW][X] = bounds[LOW][Y] = -1;
  291.     bounds[HIGH][X] = bounds[HIGH][Y] = 1;
  292.     bounds[LOW][Z] = cone->start_dist;
  293.     bounds[HIGH][Z] = 1;
  294.     /*
  295.      * Transform bounding box to world space.
  296.      */
  297.     BoundsTransform(&cone->trans.trans, bounds);
  298. }
  299.  
  300. char *
  301. ConeName()
  302. {
  303.     return coneName;
  304. }
  305.  
  306. void
  307. ConeStats(tests, hits)
  308. unsigned long *tests, *hits;
  309. {
  310.     *tests = ConeTests;
  311.     *hits = ConeHits;
  312. }
  313.  
  314. void
  315. ConeMethodRegister(meth)
  316. UserMethodType meth;
  317. {
  318.     if (iConeMethods)
  319.         iConeMethods->user = meth;
  320. }
  321.